Testing derivaties for 1D MT problem.
Especially the rx.projectFieldsDeriv
In [1]:
import SimPEG as simpeg
import simpegEM as simpegem, simpegMT as simpegmt
from SimPEG.Utils import meshTensor
import numpy as np
In [2]:
simpegmt.FieldsMT.FieldsMT_1D
Out[2]:
In [3]:
# Setup the problem
sigmaHalf = 1e-2
# Frequency
nFreq = 33
# freqs = np.logspace(3,-3,nFreq)
freqs = np.array([100])
# Make the mesh
ct = 5
air = meshTensor([(ct,25,1.3)])
# coreT0 = meshTensor([(ct,15,1.2)])
# coreT1 = np.kron(meshTensor([(coreT0[-1],15,1.3)]),np.ones((7,)))
core = np.concatenate( ( np.kron(meshTensor([(ct,15,-1.2)]),np.ones((10,))) , meshTensor([(ct,20)]) ) )
bot = meshTensor([(core[0],10,-1.3)])
x0 = -np.array([np.sum(np.concatenate((core,bot)))])
# Change to use no air
m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core))], x0=x0)
# Make the model
sigma = np.zeros(m1d.nC) + sigmaHalf
sigma[ m1d.gridCC > 0 ] = 1e-8
rxList = []
for rxType in ['z1dr','z1di']:
rxList.append(simpegmt.SurveyMT.RxMT(simpeg.mkvc(np.array([0.0]),2).T,rxType))
# Source list
srcList =[]
tD = False
if tD:
for freq in freqs:
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1DhomotD(rxList,freq))
else:
for freq in freqs:
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq))
# Make the survey
survey = simpegmt.SurveyMT.SurveyMT(srcList)
# Set the problem
problem = simpegmt.ProblemMT1D.eForm_psField(m1d)
problem.sigmaPrimary = sigma
problem.pair(survey)
# Get the fields
fields = problem.fields(sigma)
# Project the data
data = survey.projectFields(fields)
In [4]:
%debug
We need calculate this derivative. \begin{align} \underbrace{\frac{\partial P(f(u(m)),m^{fix})}{\partial f}}_{Rx} \end{align}
Use the rule \begin{align} \frac{d}{dx}\left( \frac{a(x)}{b(x)} \right) = \frac{\frac{d }{dx} a(x) b(x) - a(x)\frac{d }{dx} b(x) }{ b(x)^2 } \end{align}
In the case of the 1D MT problem the data is calculated as \begin{align} MT1Ddata = P(f(m)) &= \frac{P_{ex} f_e(src,m)}{P_{bx} f_b(src,m) \frac{1}{\mu_0}} = \frac{P_e u}{P_b f_b(u)} \\ \frac{\partial P(f(m))}{\partial u} v &= \frac{P_e}{P_b \frac{1}{mu_0} f_b(u)}v - \frac{P_e u}{\left(P_b \frac{1}{mu_0} f_b(u)\right)^2} P_b \frac{1}{mu_0} \frac{d f_b}{du} v \end{align} where u is the fields that we solve for. \begin{align} \frac{d f_b}{du} = - \frac{1}{i \omega} \nabla \end{align}
In [5]:
# Unused code &= \frac{ P_{ex} P_{bx} \frac{1}{\mu_0} \left( f_b(src,m) - f_e(src,m) \right) } { \left(P_{bx}f_b(src,m) \frac{1}{\mu_0} \right)^2 }
As matrices the formulas above can be written as \begin{align} \left[ \frac{\partial P(f(m))}{\partial u} v \right] = \left[ diag \left[ \frac{1}{\left(P_b \frac{1}{mu_0} f_b(u)\right)} \right] [P_e v] , diag[P_e u] diag \left[ \frac{1}{\left(P_b \frac{1}{mu_0} f_b(u)\right)} \right]^T diag \left[ \frac{1}{\left(P_b \frac{1}{mu_0} f_b(u)\right)} \right] \left[ P_b \frac{1}{mu_0} \frac{d f_b}{du}(v) \right] \right] \end{align}
The adjoint problem is done simliarly \begin{align} \left[ \frac{\partial P(f(m))}{\partial u} v \right]^T = [P_e]^T diag \left[ \frac{1}{\left(P_b \frac{1}{mu_0} f_b(u)\right)} \right]^T v - \left[ P_b \frac{d f_b}{du} \frac{1}{mu_0} \right]^T diag \left[ \frac{1}{\left(P_b \frac{1}{mu_0} f_b(u)\right)} \right] diag \left[ \frac{1}{\left(P_b \frac{1}{mu_0} f_b(u)\right)} \right]^T diag \left[ P_e u \right]^T v \end{align}
In [6]:
# def projectFields(self, src, mesh, u):
# '''
# Project the fields and return the
# '''
# if self.projType is 'Z1D':
# Pex = mesh.getInterpolationMat(self.locs,'Fx')
# Pbx = mesh.getInterpolationMat(self.locs,'Ex')
# ex = Pex*mkvc(u[src,'e_1d'],2)
# bx = Pbx*mkvc(u[src,'b_1d'],2)/mu_0
# f_part_complex = ex/bx
# real_or_imag = self.projComp
# f_part = getattr(f_part_complex, real_or_imag)
# return f_part
In [7]:
# Initate things for the derivs Test
src = survey.srcList[0]
rx = src.rxList[0]
v = np.random.randn(m1d.nN)
v0 = np.random.randn(m1d.nF+m1d.nE)
u0 = np.random.randn(m1d.nN)+np.random.randn(m1d.nN)*1j
f0 = problem.fieldsPair(m1d,survey)
f0[src,'e_1dSolution'] = u0
# f0[src,'b_1d'] = -1/(1j*simpegem.Utils.EMUtils.omega(src.freq))*m1d.nodalGrad*u0
In [8]:
# Run a test
def fun(u):
f = problem.fieldsPair(m1d,survey)
f[src,'e_1dSolution'] = u
return rx.projectFields(src,m1d,f), lambda t: rx.projectFieldsDeriv(src,m1d,f0,t)
simpeg.Tests.checkDerivative(fun,u0,num=5,plotIt=False)
Out[8]:
In [9]:
rx.projectFieldsDeriv(src,m1d,f0,u0)
Out[9]:
In [10]:
rx.projectFields(src,m1d,f0)
Out[10]:
In [11]:
# Test the Jvec derivative.
In [12]:
# print '%s formulation - %s' % (fdemType, comp)
CONDUCTIVITY = 0.01
m0 = np.log(np.ones(problem.mesh.nC)*CONDUCTIVITY)
# mu = np.log(np.ones(problem.mesh.nC)*MU)
if True:
m0 = m0 + np.random.randn(problem.mesh.nC)*CONDUCTIVITY*1e-1
# mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1
# prb.mu = mu
# survey = prb.survey
def fun(x):
return survey.dpred(x), lambda x: problem.Jvec(m0, x)
simpeg.Tests.checkDerivative(fun, m0, num=4, plotIt=False)
Out[12]:
In [13]:
%debug
In [14]:
### Adjoint test
We have \begin{align} Jvec =&
In [15]:
# Run a test
TOL = 1e-4
FLR = 1e-20
def projectFieldsAdjointTest(fdemType, comp):
print 'Adjoint %s formulation - %s' % (fdemType, comp)
m = np.log(np.ones(problem.mesh.nC)*0.01)
if True:
m = m + np.random.randn(problem.mesh.nC)*0.01*1e-1
u = problem.fields(m)
v = np.random.randn(1)#+np.random.randn(1)*1j
# print prb.PropMap.PropModel.nP
w = np.random.randn(m1d.nN)+np.random.randn(m1d.nN)*1j
vJw = v.dot(rx.projectFieldsDeriv(src,m1d,f0,w))
wJtv = w.dot(rx.projectFieldsDeriv(src,m1d,f0,v,adjoint=True)).real
tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR])
print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
return np.abs(vJw - wJtv) < tol
projectFieldsAdjointTest('e','projectFieldsDeriv')
Out[15]:
In [21]:
# Run a test
TOL = 1e-4
FLR = 1e-20
def getADeriv_mAdjointTest():
print 'Adjoint test e formulation - getADeriv_m'
m = np.log(np.ones(problem.mesh.nC)*0.01)
if True:
m = m + np.random.randn(problem.mesh.nC)*0.01*1e-1
u = problem.fields(m)
v = np.random.randn(m1d.nN)#+np.random.randn(1)*1j
# print prb.PropMap.PropModel.nP
w = np.random.randn(m1d.nC)#+np.random.randn(m1d.nN)*1j
vJw = v.dot(problem.getADeriv_m(freq,f0,w))
wJtv = w.dot(problem.getADeriv_m(freq,f0,v,adjoint=True))
tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR])
print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
return np.abs(vJw - wJtv) < tol
getADeriv_mAdjointTest()
In [17]:
%debug
In [18]:
# Run a test
TOL = 1e-4
FLR = 1e-20
def getRHSDeriv_mAdjointTest():
print 'Adjoint test e formulation - getRHSDeriv_m'
m = np.log(np.ones(problem.mesh.nC)*0.01)
if True:
m = m + np.random.randn(problem.mesh.nC)*0.01*1e-1
u = problem.fields(m)
v = np.random.randn(m1d.nN)#+np.random.randn(1)*1j
# print prb.PropMap.PropModel.nP
w = np.random.randn(m1d.nC)#+np.random.randn(m1d.nN)*1j
vJw = v.dot(problem.getRHSDeriv_m(freq,w))
wJtv = w.dot(problem.getRHSDeriv_m(freq,v,adjoint=True))
tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR])
print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
return np.abs(vJw - wJtv) < tol
getRHSDeriv_mAdjointTest( )
Out[18]:
In [ ]:
simpeg.mkvc(np.random.randn(survey.nD)+np.random.randn(survey.nD)*1j,2)
print survey.nD
In [19]:
TOL = 1e-4
FLR = 1e-20
def JvecAdjointTest():
print 'Adjoint e formulation - Jvec'
m = np.log(np.ones(problem.mesh.nC)*0.01)
if True:
m = m + np.random.randn(problem.mesh.nC)*0.01*1e-1
u = problem.fields(m)
v = np.random.rand(survey.nD)
# print prb.PropMap.PropModel.nP
w = np.random.rand(problem.mesh.nC)
vJw = v.dot(problem.Jvec(m, w, u))
wJtv = w.dot(problem.Jtvec(m, v, u))
tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR])
print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
return np.abs(vJw - wJtv) < tol
In [20]:
JvecAdjointTest()
Out[20]:
In [ ]: